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Performance of improved high-order filter 
schemes for turbulent flows with shocks 

By D. Kotov, H.C. Yee & B. Sjogreen 


1. Motivation and Objectives 

In simulations of multiscale compressible turbulent flows with discontinuities, the con- 
structions of numerical schemes for stable and accurate simulation of shock/turbulence 
interaction share one important ingredient - minimization of numerical dissipation while 
maintaining numerical stability. Methods commonly used for these simulations rely on 
switching between spectral or high-order compact schemes to shock-capturing schemes 
and are not always practical for multiscale shock/turbulence interactions. Also frequent 
switching between these two types of schemes can create severe numerical instability. 
One possible way of overcoming these difficulties is by using high-order nonlinear filter 
schemes (Yee & Sjogreen 2007, 2009; Sjogreen & Yee 2004). These schemes take advan- 
tage of the effectiveness of the nonlinear dissipation contained in good shock-capturing 
schemes as stabilizing mechanisms at locations where needed. In addition to the mini- 
mization of numerical dissipation while maintaining numerical stability in compressible 
turbulence with strong shocks, Yee & Sjogreen (2002); Yee (2002); Yee & Sweby (1997); 
Yee et al. (1999) discussed a general framework for the design of such schemes. 

One of the key parameters responsible for minimizing numerical dissipation in nonlin- 
ear filter schemes is the numerical dissipation control parameter n. In the original version 
of filter schemes this parameter required some tuning based on the knowledge of the flow 
physics for the particular case. The technique considered in this study made use of a k 
dependent on flow location that is a function of the local Mach number M in such a 
way that the scheme would be suitable for a wide range of flow types without additional 
tuning. The basic idea of developing a local flow sensor with improved control on nu- 
merical dissipation that is inherited in the chosen shock-capturing scheme was proposed 
in Yee & Sjogreen (2009). In particular, a local Mach sensor k(M) was proposed by 
Yee & Sjogreen without extensive numerical experiments. This is a sequel study to that 
work with comparative evaluation for the same test cases as in Johnsen et al. (2010) to 
illustrate the performance of the local Mach Sensor. The results obtained using constant 
k based on the free stream Mach number have been presented in Johnsen et al. (2010). 
This study addresses the same problems as in Johnsen et al. (2010) and compares the 
results obtained using k = const with the new results obtained by applying global or 
local k numerical dissipation control. 

2. Improved High-Order Filter Schemes 

2.1. Original High- Order Filter Scheme 
Consider the 3-D compressible Euler equations in Cartesian geometry, 
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Here the velocity vector u = (u,v,w) T , p is the density, e is the total energy, and p is 
the pressure. 

The high-order nonlinear filter scheme of Yee et al. (2007) when used to solve the fully 
coupled system (2.1), consists of three steps. 

2.1.1. Preprocessing Step 

Before the application of a high-order non-dissipative spatial base scheme, a prepro- 
cessing step is employed to improve the stability. The inviscid flux derivatives of the 
governing equation(s) are split in the following three ways, depending on the flow types 
and the desire for rigorous mathematical analysis or physical argument. 

• Entropy splitting of Olsson & Oliger (1994) and Yee et al. (2000); Yee & Sjogreen 
(2002). The resulting form is non-conservative and the derivation is based on entropy 
norm stability with boundary closure for the initial value boundary problem. 

• The system form of the Ducros et al. (2000) splitting. This is a conservative splitting 
and the derivation is based on physical arguments. 

• Tadmor entropy conservation formulation for systems Sjogreen & Yee (2009). The 
derivation is based on mathematical analysis. It is a generalization of Tadmor’s entropy 
formulation to systems and has not been fully tested on complex flows. 

2.1.2. Base Scheme Step 

A full time step is advanced using a high-order non-dissipative (or very low dissipation) 
spatially central scheme on the split form of the governing partial differential equations 
(PDEs). Summation-by-parts (SBP) boundary operator Olsson (1995); Sjogreen & Yee 
(2007) and matching order conservative high-order free stream metric evaluation for 
curvilinear grids Vinokur & Yee (2002) are used. High-order temporal discretization 
such as the third-order or fourth-order Runge-Kutta (RK3 or RK4) method is used. It 
is remarked that other temporal discretizations can be used for the base scheme step. 

2.1.3. Post-Processing (Nonlinear Filter Step) 

After the application of a non-dissipative high-order spatial base scheme on the split 
form of the governing equation(s), to further improve nonlinear stability from the non- 
dissipative spatial base scheme, the post-processing step is used to nonlinearly filter the 
solution by a dissipative portion of a high-order shock-capturing scheme with a local 
flow sensor. The flow sensor provides locations and amounts of built-in shock-capturing 
dissipation that can be further reduced or eliminated. At each grid point, a local flow 
sensor, e.g., a multi-resolution wavelet, would be employed to analyze the regularity of 
the computed flow data. Only the discontinuity locations would receive the full amount 
of shock-capturing dissipation. In smooth regions, no shock-capturing dissipation would 
be added. In regions with strong turbulence, a small fraction of the shock-capturing 
dissipation would be added to improve stability. The nonlinear dissipative portion of a 
high-resolution shock-capturing scheme can be any shock-capturing scheme. 

Let U* be the solution after the completion of the base scheme step. The final update 
of the solution after the filter step is (with the numerical fluxes in the y- and z-directions 
suppressed as well as their corresponding y- and z-directions indices on the x inviscid 
flux suppressed) 

U?ti = - |^+i/2 - H U/2^ H >+i/ 2 = Rj+HiHj+112, (2.2) 

where Rj+ 1/2 is the matrix of right eigenvectors of the Jacobian of the inviscid flux vector 
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in terms of the Roe’s average states based on U* . H* +1 / 2 and //*_-, , 2 are numerical fluxes 
in terms of the Roe’s average states based on U* . Denote the elements of the vector 
Hj. |_i /2 by hj +1 j 2 ,l = 1, 2, 5. The nonlinear portion of the filter hj +1 / 2 has the form 


h j + 1/2 - 


9 W j + l/2^ + l/2' 


(2.3) 


Here cA +1 , 2 is the wavelet flow sensor to activate the nonlinear numerical dissipation 
2 0 j _ i_i / 2 anc ^ th e original formulation for re is a positive parameter that is less than or 
equal to one. The choice of the parameter re can be different for different flow types. In the 
next section the local and the global implementation of re is discussed. These techniques 
of choosing re can be used to avoid tuning scheme parameters for a wide range of flows. 


2.2. An efficient global re for low Mach number and smooth flows 

To overcome the shortcomings of “low speed Roe scheme”, the flow speed indicator 
formula of Li & Gu (2008) was modified to obtain an improved global re denoted by re 
for (2.3) to minimize the tuning of the original re for low Mach number flows, re has the 
form: 


« = fl(M)K, 


(2.4) 


with 


fi(M) 


min 


I'M 2 v /4+(l-M 2 ) 2 
V~2 1 + M 2 



(2.5) 


Here M is the maximum Mach number of the entire computational domain at each stage 
of the time evolution. fi(M) has the same form as in Li & Gu (2008) except for the extra 
factor “M/2” added to the first argument on the right-hand-side of the original form 
f(M) in equation (18) of Li & Gu (2008). The added factor provides a similar value of 
the tuning re observed from numerical experimentation reported in aforementioned works. 
With the flow speed indicator /i(M) in front of re, the same re used for the supersonic 
shock problem can be used without any tuning for the very low speed turbulent flow 
cases. Another minor modification of the above is, 


fi(M) = max 


min 


fM 2 v/4+(l-M 2 ) 2 
V~2 1 + M 2 



( 2 . 6 ) 


where e is a small threshold value to avoid completely switching off the dissipation. A 
function which retains the majority of fi(M) but includes larger Mach number for not 
weak shocks is 


/ 2 (M) = (Q(M,2) + Q(Af,3.5))/2 


(2.7) 


or 


/ 2 (M) = max((Q(M, 2) + Q(M, 3.5))/2, e), (2.8) 


where 


Q(M, a) 


P{M/a) M < a 
1 otherwise 


(2.9) 


The polynomial 

P{x)=x 4 ( 35 - 84a; + 70a; 2 - 2Cte 3 ) (2.10) 

is monotonically increasing from P( 0) = 0 to P(l) = 1 and has the property that P'(x) 
has three continuous derivatives at x = 0 and at x = 1. 
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Figure 1 . Mach sensors: original function / of Li and Gu and its modifications /i and /2 

Below supersonic speeds a simple and efficient global 7t can be obtained according to 
the maximum Mach number of the entire flow field and the value is determined by fi(M) 
or f 2 {M) for non-zero w*- +1 / 2 - It is noted that if the original /(M) were used instead of 
/i(M) or f 2 (M) in Eq.(2.4), the amount of nonlinear filter dissipation could be too large 
for very low speed turbulent flows (for the same fixed n). See Fig. 1 for details. 

2.3. Local flow sensor for a wide spectrum of flow speed and shock strength 

At each time step and grid point, the aforementioned global ~k is not sufficient to reduce 
the amount of numerical dissipation where needed for flows that contains a variety of 
flow features. A more appropriate approach is to obtain a “local k” that is determined 
similar to K at each grid point. If known, a dominating shock jump variable should be 
used for shock detections. In other words, the filter numerical flux indicated in Eq.(2.3) 
is replaced by: 

hj + l/2 = 2 Wj+l/2 u, j+\/2 ( t ) j+l/2\- (2-11) 

In this formula all the variables are supposed to be evaluated at the Roe-average 
state. However, the results presented in this paper are obtained with simple averaging 
for local K l j + i / 2 , i- 0 . , Kj+ 1/2 = (/Cj + Kj+ 1)/2 for efficiency. However, this approach was 
not sufficient for the case with isotropic turbulence, where we use Kj+ 1/2 = max(/tj, Kj+i) 
over all 3 directions. 

In the case of unknown physics and without experimental data or theory for compari- 
son, K , j + - l /2 has to depend on the local Mach number in low speed or smooth flow regions, 
on local shock strength in shock regions and on turbulent fluctuations in vortical regions 
in order to minimize the tuning of parameters. According to the local flow type, for each 
non-zero wavelet indicator Wj +1 / 2 , k^ + 1 , 2 should provide the aforementioned amount 
(between (0,1)) to be filtered by the shock-capturing dissipation c ■ For problems 
containing turbulence and strong shocks, the shock strength should come into play. One 
measure of the shock strength can be based on the numerical Schlieren formula Hadjadj 
& Kudryavtsev (2005) for the chosen variables that exhibit the strongest shock strength. 
In the vicinity of turbulent fluctuation locations, k*. + 1 , 2 will be kept to the same order 
as in the nearly incompressible case, except in the vicinity of high shear and shocklets. 


5 


Improved high-order filter schemes 

Due to the fact that k works well for local Mach number below 0.4, n only needs to 
be modified in regions that are above 0.4. In other words, the final value of k*- + 1 / 2 is 
determined by the previous local k, if the local Mach number is below 0.4. If the local 
Mach number is above 0.4, at discontinuities detected by the non-zero wavelet indicator, 
w j+i/2’ K j+ 1/2 is determined by the shock strength (normalized between [0,1]) based 
on the Schlieren formula near discontinuities. Again, if known, dominating shock jump 
variables should be used for shock detection. 

At locations with turbulence, determined by the turbulent sensor (e.g., ur. +1 / 2 obtained 
from employing wavelets with higher order vanishing moments), Hj +1 / 2 is kept to the 
same order as in the nearly incompressible case, except in the vicinity of high shear and 
shocklet locations, where a slightly larger k*. + 1 , 2 would be used. Wavelets for detecting 
turbulent flow can be (a) Wavelets with higher order vanishing moments, and (b) Wavelet 
based Coherent Vortex Extraction (CVE) of Farge et al. (2001). Since most flow sensors 
are not optimal for all flow types, other sensors were discussed in Yee & Sjogreen (2009). 
For example, in a multiscale flow, the best fit local flow sensor would require the switching 
among (or the combination of) wavelets, Ducros et al. (Ducros et al. 2000) and gradient 
sensors (Ducros et al. 2000) from region to region of the flow field. 


3. Test cases 

In this section we illustrate the performance of the high-order filter schemes with 
different dissipation control parameters on a set of benchmark solutions. The first three 
problems are inviscid and the last one is viscous. For each test case the results are 
obtained on a coarse grid and compared with the reference solution, which can be either 
an analytical approximation or a solution by a standard method using a grid convergence 
process. Also we compare the performance of the dissipation control using local k with 
the one using the global k. Note that the results obtained by filter methods presented in 
Johnsen et al. (2010) are using constant k, which depending on the choice of the constant, 
may obtain similar results to the case when using global k technique. Also, all test cases 
include the results obtained by the filter scheme with constant dissipation control for 
which only values of re = 1 and k = 0.4 have been picked in order not to complicate the 
plots. 

For comparison with the the standard schemes WENO of 7 th order is used, designated 
as “WEN07”. All the filter schemes considered here imply central schemes of 8 th order 
with Ducros splitting filtered with the dissipative portion of WEN07 and designated as 
“WEN07fi+split” . 


3.1. Shu- Osher problem 

The first test case we consider is the Shu-Osher problem (Shu & Osher 1989). This is a 
one-dimensional idealization of shock-turbulence interaction in which a shock propagates 
into a perturbed density field. The goal of this problem is to test the capability to 
accurately capture a shock wave, its interaction with an unsteady density field, and the 
waves propagating downstream of the shock. The one-dimensional Euler equations with 
7 = 1.4 are solved on the domain x € [—5, 5] with Ax = 0.05 and initial conditions 


( (3.857143, 2.629369, 10.33333) x < -4, 
\ (l + 0.2sin(5x),0,l) x > 4. 


(3.1) 


This problem corresponds to a M = 3 shock moving into a field with a small density 
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Figure 2. Density profiles for the Shu-Osher problem with Ax = 0.05 and t = 1.8. Compari- 
son of WEN07fi+split with different control of dissipation amount. Reference solution is with 
Ax = 6.25 x 10 -4 . 


(or entropy) disturbance. The solution is compared to a reference solution with Ax = 
6.25 x 1CT 4 . 

The density and velocity profiles at time t = 1.8 are shown on Figs 2 and 3. The 
interaction between the shock and the entropy disturbance generates both acoustic and 
entropy waves downstream of the shock. The acoustic waves are strong enough to steepen 
into weak shock waves. At the given time, the shock location is x s ~ 2.39, the location 
of the contact discontinuity at the leading entropy wave is x c ps 0.69, and the location of 
the leading acoustic wave is x a ss 2.75. 

For this test case the results obtained by the scheme using k = 0.4 contain spurious 
oscillations because of too small dissipation. On the contrary, the results obtained by 
regular WEN07 are too dissipative. The results of all the other schemes are close to 
each other and in better agreement with the reference solution. The scheme with local 
k contains a little more dissipation but fewer oscillations than the schemes with global 
k and k = 1. Note that when k = 1, the full amount of shock-capturing numerical 
dissipation is being employed. 

3.2. Shock-vorticity/ entropy wave interaction 

A generalization of the Shu-Osher problem is the two-dimensional interaction of a vor- 
ticity/entropy wave with a normal shock (Mahesh 2000). The two-dimensional Euler 
equations are solved with 7 = 1.4 on the domain x € [0, 47 r],y £ [— 7 r, 7 r) with Ax = 
7r/50,A y = 7t/ 16. Periodic boundaries are used in the y-direction; x = 0 is a super- 
sonic inflow and x = 47 r is a subsonic outflow. Different techniques are employed to 
avoid acoustic reflections from the outflow, including an extension of the domain with a 
sponge region. First, a one-dimensional base solution corresponding to a M = 1.5 shock 
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Reference 

O WEN07 

WEN07fi+split, k = 1 

WEN07fi+split, k = 0.4 

-WEN07fi+split, global k 
WEN07fi+split, local k 



Figure 3. Zoom of density (left) and velocity (right) profiles for the Shu-Osher problem with 
Ax = 0.05 and t = 1.8. Comparison of WEN07fi+split with different control of dissipation 
amount. Reference solution is with Ax = 6.25 x 10 -4 . 



Figure 4. Instantaneous vorticity contours for the shock- vorticity/entropy wave interaction at 
t = 25 obtained by WEN07fi-|-split. Left: global dissipation control, right: local dissipation 
control. 


is defined as 

/- - f (pl,ul,Pl) = (1) 1-5, 0.714286), x<3ir/2 1 

[P,u x ,p) | ^ Pr Ur Pr ^ _ (i + 0.2sin(5x),0, 1), x > 3tt/2. 


(3.2) 


A combined vorticity/entropy wave is superposed onto the base flow. The initial data 
then becomes 


P = P + PlA s cos(k x x + k y y), 

u x = u x + ulA v sini/i cos(k x x + k y y), 

u y = — ulA v cos ip cos(k x x + k v y), 

p = p. 

The conditions at the inflow boundary x = 0 are: 

p = Pl+ PlA c cos [k v y - k x u L t), 
u x =u L + u l A v sin xp cos(k y y - k x u L t), 
u y = —ulA v cos tp cos (k y y — k x UL,t ), 

P = PL- 


(3.3) 


(3.4) 


For the present study, k x = 1/ tan^, %p = 45°, k y = 1, A e = A v = 0.025. 

Figure 4 shows instantaneous vorticity contours obtained by schemes with global and 
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Figure 5. Kinetic energy profiles (left) and mean-square vorticity (right) for the shock- vortic- 
ity/entropy wave interaction averaged in span and time period, obtained by WEN07fi+split 
with different control of dissipation amount. The analytical solution is the linear analysis of 
Mahesh (2000). 


local dissipation control. The results obtained with local k have less oscillations especially 
in the vicinity of the shock. Figure 5 shows kinetic energy profiles and mean-square vor- 
ticity profiles averaged in span and time period obtained by different methods. The 
vorticity profiles are compared with the linear analysis of Mahesh (2000). All methods 
show oscillatory behavior near the shock. The case with n = 0.4 has the highest oscilla- 
tions. For this case only methods with local k and with k = 1 obtain vorticity close to 
analytical level. The method with global k overestimates vorticity and the method with 
k = 0.4 underestimates it. The reason why global n control shows worse results for this 
case might be connected with using the Mach curve 2.7 for the whole range of M, while 
local k algorithm implies using it only for M < 0.4 and computing shock strength using 
the Schlieren formula. 


3.3. Taylor-Green vortex 

From a well-resolved initial condition, the inviscid Taylor-Green vortex (Taylor & Green 
1937) begins stretching and producing ever smaller scales. It thus constitutes a non- 
regularized problem with no lower bound on the length scale and is solved with no regu- 
larization other than that provided by the numerical method. The goal of this problem is 
to provide a test of the stability of the methods for severely under-resolved motions, as 
well as a measure of the preservation of kinetic energy and the growth of enstrophy. The 
three-dimensional Euler equations are solved with gas constant 7 = 5/3. The domain is 
a cube with edge length 27r and the grid size is 64 3 . Boundary conditions are periodic in 
all directions. The initial conditions are 

P = !, 

u x = sin x cos y cos z, 

Vty —— cos x sin y cos z, (3-5) 

u z = 0, 

p = 100 + ([cos(2z) + 2][cos(2x) + cos(2y)] — 2) / 16. 

where the mean pressure is sufficiently high to make the problem essentially incompress- 
ible. In the incompressible problem the kinetic energy should remain constant while the 
enstrophy grows rapidly. Figure 6 shows the results obtained by different methods. The 
left part of Fig. 6 represents the kinetic energy time evolution averaged by the computa- 
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Figure 6. Kinetic energy (left) and enstrophy (right) for Taylor-Green vortex problem obtained 
by WEN07fi+split with different control of dissipation amount. The enstrophy is also compared 
with the results from Johnsen et al. (2010) obtained by “Hybrid” code. 


tional domain. The kinetic energy is conserved only in the case of using local n. The other 
methods are not able to conserve kinetic energy, though the method with global k has 
much smaller loss of energy than methods with k = const. The right part of Fig. 6 rep- 
resents enstrophy evolution obtained by filter schemes and “Hybrid” code from Johnsen 
et al. (2010). Inside the range t € [0; 6.5], both methods with global and local dissipation 
control exhibit high accuracy compare with results published in Johnsen et al. (2010). 
Both methods with k = 1 and n = 0.4 underestimate the growth of the enstrophy. 


3.4. Compressible isotropic turbulence 

The final test case is that of decaying compressible isotropic turbulence with eddy shock- 
lets (Lee et al. 1991). Given a sufficiently high turbulent Mach number M t , weak shock 
waves (eddy shocklets) develop spontaneously from the turbulent motions. The goal of 
this problem is to test the ability of the methods to handle shocklets when locations 
are not known a priori, as well as the accuracy for broadband motions in the presence 
of shocks. The three-dimensional Navier-Stokes equations with 7 = 1.4 are solved on 
the cube with the edge size 27r, grid size 64 3 and periodic boundary conditions in all 
directions. The physical viscosity is assumed to follow a power-law of the type: 

p/Pref = ( T/T ref f /A (3.6) 


The initial condition consists of a random solenoidal velocity field Ui t 0 that satisfies 


E(k) ~ fc 4 exp(— 2(fc/fco) 2 ), ^rms,o = E(k)dk (3.7) 


The brackets here denote averaging over entire computational domain. For this study, 
in order to be consistent with Johnsen et al. (2010), we choose u rms ,o = 1 and = 4. 
The density and pressure fields are initially constant with initial turbulent Mach number 
M t p = 0.6 and Taylor-scale Reynolds Re\p = 100. These parameters are defined as 
follows: 


M t = 


<s/(UiUi) 

(c) ’ 


Re\ = 


(p)< 


(p) 


(3.8) 
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Figure 7. Temporal evolution of the variance of different quantities for the isotropic turbulence 
problem on 64 3 grid. Top row: left - kinetic energy, right - enstrophy. Bottom row: left - tem- 
perature variance, right - dilatation, = diUi. The reference is the solution from Johnsen et al. 
(2010) on a 256 3 grid spectrally filtered to a 64 3 grid. 


where 


'U'rms 




(UiUi) 

3 


A = 


I (ug) 

{(d x u x ) 2 ) ' 


(3.9) 


Time scale is r = Xo/u rms fi and the final time is t/r = 4. 

The temporal evolution of the mean-square velocity and vorticity, the temperature 
variance and dilatation are shown in Fig. 7. The reference solution is the solution from 
Johnsen et al. (2010) obtained by DNS on 256 3 grid spectrally filtered to a 64 3 grid. 
All four methods with dissipation control are much closer to the reference solution than 
standard WEN07, and all these four methods agree quite well with DNS results for 
mean-square velocity and temperature variance. However, only the method with local k 
agrees well with DNS results for vorticity (right-top plot Fig. 7). Also this method is the 
closest one to DNS dilatation results (right bottom plot on Fig. 7). 

Instantaneous profiles of dilatation, density, temperature variance and vorticity through 
an eddy shocklet are shown in Fig. 8. The shocklet is located at x « 2.8, i.e. , at the point 
of large negative dilatation. All methods show reasonable results for such coarse grid as 
64 3 and as expected none of them could resolve the shocklet. The biggest discrepancy 
from DNS results is obtained by regular WEN07 and all the other four methods obtained 
similar results. The method with local k obtains results slightly closer to the reference 
then other methods. 
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Figure 8. Instantaneous profiles through a shocklet (at x « 2.8 ) for the isotropic turbulence 
problem on 64 3 grid at t/r = 4. Top row: left - kinetic energy, right - enstrophy. Bottom row: left 
- temperature variance, right - dilatation, 9i = diUi. The reference is the solution from Johnsen 
et al. (2010) on a 256 3 grid spectrally filtered to a 64 3 grid. 

4. Conclusions 

The performance of the filter scheme with improved dissipation control k has been 
demonstrated for different flow types. The scheme with local k is shown to obtain more 
accurate results than its counterparts with global or constant k. At the same time no 
additional tuning is needed to achieve high accuracy of the method when using the local 
k technique. However, further improvement of the method might be needed for even more 
complex and/or extreme flows. 
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